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Abstract 

In this paper we present a modeling approach to bridge the atomistic with macroscopic 
scales in crystalline materials. The methodology combines identification and modeling of the 
controlling unit processes at microscopic level with the direct atomistic determination of fun- 
damental material properties. These properties are computed using a many body Force Field 
derived from ab initio quantum-mechanical calculations. This approach is exercised to de- 
scribe the mechanical response of high-purity Tantalum single crystals, including the effect of 
temperature and strain-rate on the hardening rate. The resulting atomistically informed model 
is found to capture salient features of the behavior of these crystals such as: the dependence of 
the initial yield point on temperature and strain rate; the presence of a marked stage I of easy 
glide, specially at low temperatures and high strain rates; the sharp onset of stage II hardening 
and its tendency to shift towards lower strains, and eventually disappear, as the temperature in- 
creases or the strain rate decreases; the parabolic stage II hardening at low strain rates or high 
temperatures; the stage II softening at high strain rates or low temperatures; the trend towards 
saturation at high strains; the temperature and strain-rate dependence of the saturation stress; 
and the orientation dependence of the hardening rate. 
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1 Introduction 



This paper is concerned with the development of a multiscale modeling approach for advanced 
materials such as high-purity bcc single crystals. The present approach is aligned with the current 
divide and conquer paradigm in micromechanics (see, e. g., [[I], 0, |3|, §]. This paradigm 
first identifies and models the controlling unit process at microscopic scale. Then, the energetics 
and dynamics of these mechanisms are quantified by means of atomistic modeling. Finally, the 
macroscopic driving force is correlated to macroscopic response via microscopic modeling. This 
last step involves two stages, localization of the macroscopic driving force into unit-process driving 
forces and averaging of the contribution of each unit process into the macroscopic response. 

We show in this article that the meticulous application of this paradigm renders truly predictive 
models of the mechanical behavior of complex systems. In particular we predict the hardening of 
Ta single crystal and its dependency for a wide range of temperatures and strain rates. The feat 
of this approach is that predictions from these atomistically informed models recover most of the 
macroscopic characteristic features of the available experimental data, without a priori knowledge 
of such experimental tests. This approach then provides a procedure to forecast the mechanical 
behavior of material in extreme conditions where experimental data is simply not available or very 
difficult to collect. 

A crucial step in this approach is the appropriate selection and modeling of the unit processes. 
These models supply the link between the atomic and meso scale by identifying and correlating the 
relevant material properties, susceptible of atomistic determination such as energy formation for 
defects, with the corresponding driving forces. In this case, we specifically consider the following 
unit processes: double-kink formation and thermally activated motion of kinks; the close-range 
interactions between primary and forest dislocation, leading to the formation of jogs; the percola- 
tion motion of dislocations through a random array of forest dislocations introducing short-range 
obstacles of different strengths; dislocation multiplication due to breeding by double cross-slip; 
and dislocation pair-annihilation. 

A set of material parameters is then obtained from the modeling and identification stage, which 
is required to quantify the contribution of each of the unit processes. We compute these mate- 
rials properties using a combination of ab-initio quantum mechanics (QM) and Force Field (FF) 
calculations. QM describes the atomic interactions from first principles, i.e. with no input from 
experiments; unfortunatelly QM methods are computationally intensive and restricted to small 
systems, making QM calculations impractical to study most of the materials properties govern- 
ing plasticity. Force Fields give the total energy of a system as a potential energy function of the 
atomic positions and with Molecular Dynamics (MD) allows the simulation of systems containing 
millions of atoms. We used ab-initio quantum mechanical calculations (equations of state of vari- 
ous crystalline phases, elastic constants, energetics of defects, etc.) to develop a many body Force 
Field (FF) (named qEAM FF) for Tantalum. We use the qEAM FF with MD to calculate the core 
energy of the l/2a < 111 > screw dislocation, that of the edge dislocation with Burgers vector 
b — l/2a < 111 > in (110) planes. We have also calculated the formation energies and nucleation 
lengths of the kinks in6 = l/2a<lll> screw dislocations. 

The organization of the paper follows the sequential stages of the proposed approach. First, we 
provide a brief description of each of the unit processes including the governing final equations. 
We then identify and compute by atomistic means the corresponding material properties. Finally, 
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(a) (b) 
Figure 1 : Schematic of the double-kink mechanism. 

we compare the predictions against experimental data. 

2 Unit Processes 

Plastic deformation in metallic systems is the macroscopic manifestation of dislocation activity. 
The resistance to the dislocation motion, therefore, engenders the hardening properties observed 
in this type of materials. It is then the complex interplay of microscopic mechanisms control- 
ling dislocation mobility, dislocation interaction and dislocation evolution which confers the 
macroscopic constitutive properties. In the present approach, these controlling processes are con- 
sidered to be orthogonal in the sense that are weakly coupled with each other. The interaction 
among them is only established through the uniqueness of the macroscopic driving force which 
are shared, via the localization process, by all the unit processes. 

In this section, we introduce the set of controlling unit processes which have been identified for 
describing the mechanical response of high-purity BCC single crystals, in particular for Tantalum. 
We also provide the final expression resulting from the the modeling of each of these processes. A 
detailed description of the model, including comparison with experimental data is given in [0] . 

2.1 Dislocation Mobility: Double-Kink Formation and Thermally Activated 
Motion of Kinks 

We consider the thermally activated motion of dislocations within an obstacle-free slip plane. Un- 
der these conditions, the motion of the dislocations is driven by an applied resolved shear stress r 
and is hindered by the lattice resistance, which is weak enough that it may be overcome by thermal 
activation. The lattice resistance is presumed to be well-described by a Peierls energy function, 
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Figure 2: Temperature dependence of the effective Peierls stress for various strain rates. Note that 
the typical order of magnitude of 7g mk = 10 6 

which assigns an energy per unit length to dislocation segments as a function of their position on 
the slip plane. 

In bcc crystals, the core of screw dislocation segments relaxes into low-energy non-planar 
configurations [§, |9|, [K| [n], [12|, |5|, [13], [T4Q . This introduces deep valleys into the Peierls energy 



function aligned with the Burgers vector directions and possessing the periodicity of the lattice. At 
low temperatures, the dislocations tend to adopt low-energy configurations and, consequently, the 
dislocation population predominantly consists of long screw segments. In order to move a screw 
segment normal to itself, the dislocation core must first be constricted, which requires a substantial 
supply of energy. Thus, the energy barrier for the motion of screw segments, and the attendant 
Peierls stress, may be expected to be large, and the energy barrier for the motion of edge segments 



to be comparatively smaller. For instance, Duesbery and Xu QT5Q have calculated the Peierls stress 
for a rigid screw dislocation in Mo to be 0.022/i, where fi is the (111) shear modulus, whereas 
the corresponding Peierls stress for a rigid edge dislocation is 0.006/i, or about one fourth of the 
screw value. This suggests that the rate-limiting mechanism for dislocation motion is the thermally 
activated motion of kinks along screw segments ([|16|, |T7[ |18]]). 

At sufficiently high temperatures and under the application of a resolved shear stress r > 
0, a double-kink may be nucleated with the assistance of thermal activation (e. g., [ |T9| , |20| , |5p, 



and the subsequent motion of the kinks causes the screw segment to effectively move forward, 
Fig. |T]. Under this conditions the following expression for the effective temperature and strain-rate 
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Figure 3: Bow-out mechanism for a dislocation segment bypassing an obstacle pair 



dependent Peierls r p is obtained 



where the effective Peierls stress is given by 



jjj kink 



and the reference strain is defined as 



7 kmk = 2bpl P u D (3) 

In the preceeding equations, b is the Burgers vector, p is the dislocation density, (3 = l/k B T, k B is 
Boltzmann's constant, T is the absolute temperature, and vp> is the attempt frequency which may be 
identified with the Debye frequency to a first approximation. Also, lp is the distance between two 
consecutive Peierls valleys. For bcc crystals, lp = y^2/3a if the slip plane is {110}, lp = \J~2a, if 
the slip plane is {112}, and l P = ^8/3a if the slip plane is {123}, where a is the cubic lattice size 
Finally, E kmk is the energy of formation of a kink-pair and L kmk is the length of an incipient 
double kink. The formation energy E kmk and the length L kmk , which cannot be reliably estimated 
from elasticity since the energy is composed mostly of core region, can, however, be accurately 
computed by recourse to atomistic models as shown in section |3|. Modeling of this first unit process 
renders the first 2 material properties amenable of atomistic calculations. 

In Figure ^ the dependence of the effective Peierls stress on temperature and rate of deforma- 
tion is illustrated. The Peierls stress decreases ostensibly linearly up to a critical temperature T c , 
beyond which it tends to zero. These trends are in agreement with the experimental observations 



5 



After intersection 



Unfavorable Ju nctlon / \ 




Before intersection 



Favorable Junction 



Energy increase due 
tojog formation 



Figure 4: Schematic of energy variation as a function of a reaction coordinate during dislocation 
intersection and crossing. 

of Wasserbach and Lachenmann and Schultz [^3]]. The critical temperature T c increases with 
the strain rate. In particular, in this model the effect of increasing (decreasing) the strain rate has 
an analogous effect to decreasing (increasing) the temperature, and vice- versa, as noted by Tang et 
al. Q24|]. In the regime of very high strain-rates (7 > 10 5 s _1 ), effects such as electron and phonon 



drag become important and control the velocity of dislocations Q25|, [261] . 



2.2 Dislocation Interactions: Obctacle-Pair Strength and Obstacle Strength 

In the forest-dislocation theory of hardening, the motion of dislocations, which are the agents of 
plastic deformation in crystals, is impeded by secondary -or 'forest'- dislocations crossing the 
slip plane. As the moving and forest dislocations intersect, they form jogs or junctions of varying 
strengths [|27|, ^9], ^J], [I2|, fJl|, which, provided the junction is sufficiently short, may 



be idealized as point obstacles. Moving dislocations are pinned down by the forest dislocations 
and require a certain elevation of the applied resolved shear stress in order to bow out and bypass 
the pinning obstacles. For the case of infinitely strong obstacles, the resistance of the forest is 
provided by the strength of the obstacle pairs. This obstacle pair strength is subsequently reduced 
by considering that point obstacles composing the pair can only provide a finite strength. The 
processes imparting the pair-obstacle strength and obstacle strength are described next 

2.2.1 Obstacle-Pair Strength 

We begin by treating the case of infinitely strong obstacles. In this case, pairs of obstacles pin down 
dislocation segments, which require a certain threshold resolved shear stress s in order to overcome 
the obstacle pair. The lowest-energy configuration of unstressed dislocation segments spanning an 
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Figure 5: Schematic of jog formation during dislocation intersection 



obstacle pair is a step of the form shown as the thin line in Fig. Q Under these conditions, the bow- 
out mechanism by which a dislocation segment bypasses an obstacle pair may be expected to result 
in the configuration shown in Fig. |3] (bold line). If the edge-segment length is l e , a displacement 
da e of the dislocation requires a supply of energy equal to 2U screw da e + brp dge l e da e in order to 



overcome the Peierls resistance rp dsc and to extend the screw segments. The corresponding energy 



release is brl e da e . Similar contributions result from a displacement da s of the screw-segment of 
length l s . Retaining dominant terms the obstacle-pair strength is 



screw | 
' P 



2Jjedge 



(4) 



The obstacle-pair strength can be therefore estimated by quantifying Tp, l s and <7 cdgc . An expres- 
sion for the Peierls stress r P is given in Eq. ([]]), the distance between obstacles along the screw 
direction l s is estimated by statistics assuming a random obstacle distribution and the core energy 
per unit length in the edge direction U edge is obtained by atomistic calculations presented in the 
following sections. 



2.2.2 Obstacle Strength 

In this section we proceed to estimate the obstacle strengths which reduces the obstacle- 
pair strength described in the previuos section. The interaction between primary and sec- 
ondary dislocations may result in a variety of reaction products, including jogs and junctions 
PTj , ^, ^9p. Experimental estimates of junction strengths have been given by 

Franciosi and Zaoui [ [351 ] f° r m e twelve slip systems belonging to the family of {111} planes and 
[110] directions in fee crystals, and by Franciosi [BH] for the twenty-four systems of types {211} 



[111] and {110} [111] in bec crystals. The strength of some of these interactions has recently been 
computed using atomistic and continuum models $2Jj [28| , |], ^9|. Tang et al. have numerically 



estimated the average strength of dislocation junctions for Nb and Ta crystals 

For purposes of the present theory, we specifically concern ourselves with short-range interac- 
tions between dislocations which can be idealized as point defects. For simplicity, we consider the 
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Figure 6: Schematic of a dislocation line overcoming a junction 

case in which each intersecting dislocation acquires a jog. The energy of a pair of crossing disloca- 
tions is schematically shown in Fig. |] as a function of some convenient reaction coordinate, such 
as the distance between the dislocations. The interaction may be repulsive, resulting in an energy 
barrier, or attractive, resulting in a binding energy, Fig. |]. In the spirit of an equilibrium theory, 
here we consider only the final reaction product, corresponding to a pair of jogged dislocations 
at infinite distance from each other, and neglect the intermediate states along the reaction path. 
In addition, we deduce the strength of the obstacles directly from the energy supply required to 
attain the final state, i. e. the jog-formation energy. Despite the sweeping nature of these assump- 
tions, the predicted saturation strengths in multiple slip are in good agreement with experiment (cf 
Section 0), which lends some empirical support to the theory. 

We estimate the jog formation energy as follows. Based on energy and mobility considerations 
already discussed, we may expect the preponderance of forest dislocations to be of screw character, 
and the mobile dislocation segments to be predominantly of edge character. We therefore restrict 
our analysis to intersections between screw and edge segments. The geometry of the crossing 
process is schematically shown in Fig. |5|. Each dislocation acquires a jog equal to the Burgers 
vector of the remaining dislocation. The energy expended in the formation of the jogs may be 
estimated as 



E jogs \ bU scrcw [l-r cos 9 a/3 ] if b Q = 

al3 ~ |fr[/scrcw \2 r - cos(^ a ) - rcosO af) ] otherwise ° 

where r = u ed g e /u screw [ s the ratio of screw to edge dislocation line energies. This ratio is 
computed by atomistic calculations presented in the next section, renders a value of r = 1.77 for 
Ta. The resulting jog formation energies for the complete collection of pairs of {211} and {110} 
dislocations are tabulated in Table |l| 



A derivation entirely analogous to that leading to Eq. (jl]) yields the following expression for 
the strength of an obstacle in the slip system a produced by a forest segment in the system j3 

PKi \i% J 
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A2 A2'A3 A3'A6 A6'B2 B2"B4 B4'B5 B5'C1 C1'C3 C3"C5 C5"D1 D1"D4 D4"D6 D6" 
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A6 

A6' 

B2 

B2" 

B4 

B4' 

B5 

B5' 

CI 

CI' 

C3 

C3" 

C5 

C5" 

Dl 

Dl" 

D4 

D4" 

D6 

D6" 



— 1.0 1.0 1.0 1.01.0 1.5 1.5 
1.0 — 1.0 1.0 1.01.0 3.23.2 
1.0 1.0 — 1.0 1.01.0 2.42.4 
1.0 1.0 1.0 — 1.01.0 1.8 1.8 
1.0 1.0 1.0 1.0 — 1.0 2.42.4 
1.0 1.0 1.0 1.0 1.0 — 1.8 1.8 
1.5 1.5 1.5 1.5 1.5 1.5 — 1.0 
3.2 3.2 3.2 3.2 3.2 3.2 1.0 — 
2.4 2.4 2.4 2.4 2.4 2.4 1.01.0 
1.8 1.8 1.8 1.8 1.8 1.8 1.01.0 

2.4 2.4 2.4 2.4 2.4 2.4 1.01.0 
1.8 1.8 1.8 1.8 1.8 1.8 1.01.0 
1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 
1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 

1.5 1.5 1.5 1.5 1.5 1.5 2.42.4 
3.2 3.2 3.2 3.2 3.2 3.2 1.8 1.8 
2.4 2.4 2.4 2.4 2.4 2.4 1.5 1.5 
1.8 1.8 1.8 1.8 1.8 1.8 3.23.2 
1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 
1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 

2.4 2.4 2.4 2.4 2.4 2.4 1.5 1.5 
1.8 1.8 1.8 1.8 1.8 1.8 3.23.2 

1.5 1.5 1.5 1.5 1.5 1.5 2.42.4 
3.2 3.2 3.2 3.2 3.2 3.2 1.8 1.8 



1.5 1.5 1.51.5 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 

3.23.2 3.23.2 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 

2.42.4 2.42.4 1.5 1.5 1.5 1.5 1.5 1.5 2.4 2.4 2.4 2.4 2.4 2.4 

1.8 1.8 1.8 1.8 3.23.2 3.23.2 3.23.2 1.8 1.8 1.8 1.8 1.8 1.8 

2.42.4 2.42.4 2.4 2.4 2.4 2.4 2.4 2.4 1.5 1.5 1.5 1.5 1.5 1.5 

1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 3.23.2 3.23.2 3.23.2 

1.01.0 1.01.0 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 

1.01.0 1.01.0 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 

— 1.0 1.01.0 2.4 2.4 2.4 2.4 2.4 2.4 1.5 1.5 1.5 1.5 1.5 1.5 

1.0 — 1.01.0 1.8 1.8 1.8 1.8 1.8 1.8 3.2 3.2 3.23.2 3.23.2 

1.01.0 — 1.0 1.5 1.5 1.5 1.5 1.5 1.5 2.4 2.4 2.4 2.4 2.4 2.4 

1.01.0 1.0 — 3.23.2 3.23.2 3.23.2 1.8 1.8 1.8 1.8 1.8 1.8 

1.8 1.8 1.8 1.8 — 1.0 1.0 1.0 1.0 1.0 3.23.2 3.23.2 3.23.2 

1.8 1.8 1.8 1.8 1.0 — 1.0 1.0 1.0 1.0 3.2 3.2 3.2 3.2 3.23.2 

2.42.4 2.42.4 1.01.0 — 1.0 1.01.0 2.4 2.4 2.4 2.4 2.4 2.4 

1.8 1.8 1.8 1.8 1.0 1.0 1.0 — 1.0 1.0 1.8 1.8 1.8 1.8 1.8 1.8 

1.5 1.5 1.51.5 1.01.0 1.0 1.0 — 1.0 2.4 2.4 2.4 2.4 2.4 2.4 

3.23.2 3.23.2 1.0 1.0 1.0 1.0 1.0 — 1.8 1.8 1.8 1.8 1.8 1.8 

1.8 1.8 1.8 1.8 3.2 3.2 3.2 3.2 3.23.2 — 1.0 1.0 1.0 1.0 1.0 

1.8 1.8 1.8 1.8 3.2 3.2 3.2 3.2 3.23.2 1.0 — 1.0 1.0 1.0 1.0 

1.5 1.5 1.51.5 2.4 2.4 2.4 2.4 2.4 2.4 1.0 1.0 — 1.0 1.01.0 

3.23.2 3.23.2 1.8 1.8 1.8 1.8 1.8 1.8 1.0 1.0 1.0 — 1.0 1.0 

2.42.4 2.42.4 2.4 2.4 2.4 2.4 2.4 2.4 1.0 1.0 1.0 1.0 — 1.0 

1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.0 1.0 1.0 1.0 1.0 — 



Table 1: Normalized jog-formation energies resulting from crossings of bcc dislocations. 
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where the strength at zero temperature is given by 



and the reference strain rate by 

7o a = 2p a bl a u D (8) 

The lengths l a and L junct describe the geometry of the junction as illustrated in Fig. @. These 
values, which have been estimated to be of the order of few b in the present case, can also be 
obtained by atomistic models. 



2.3 Dislocation Evolution: Multiplication and Attrition 

The density of forest obstacles depends directly on the dislocation densities in all slip systems 
of the crystal. Therefore, in order to close the model we require a equation of evolution for the 
dislocation densities. Processes resulting in changes in dislocation density include production by 
fixed sources, such as Frank-Read sources, breeding by double cross slip and pair annihilation 
(see [0 for a review; see also © © 13 0> HtD- Although the operation of fixed Frank- 
Read sources is quickly eclipsed by production due to cross slip at finite temperatures, it is an 
important mechanisms at low temperatures. The double cross-slip, fixed Frank-Read sources and 
pair annihilation mechanisms are next considered in turn. 



2.3.1 Dislocation Multiplication: Fixed Frank- Reed and Breeding by Cross Glide 

The rate of dislocation multiplication in a given slip system a produced by fixed Frank-Reed 
sources and by breeding by cross glide is written as 

bp a = A 0V W (9) 

where Ao is a constant associated with the fixed Frank-Read production; this parameter is rather 
topological than material dependent. 



2.3.2 Attrition: Pair Annihilation 

The rate of dislocation attrition due to pair annihilation may finally be estimated as: 

bp a = -K P a ^ a (10) 

where k is the effective annihilation distance. This is the maximum distance at which two screw 
segments with opposite direction and forced to move with a velocity v = j/bp will annihilate. 
This distance can be estimated by simply equating the time required for trapping and escaping. 
Trapping is governed by the elastic interaction forces (attraction) while escaping by the applied 
strain rate. Then, 

1 1 1 

_ = _+ 7=f^=t (ID 

K K-c K (A + y/A 2 + 1) 
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where 



A = e-^ j ° s /3£ jog 7i°77 a (12) 

is a factor depending on the strain rate and temperature, 

7^ = 2bpl P u D (13) 

is a reference slip-strain rate and k c is the cut-off value corresponding the effective screening dis- 
tance. It follows that the critical pair-annihilation distance k decreases with increasing strain rate 
and decreasing temperature. Thus, at high strain rates the dislocation velocities are high and the 
probability of being captured by another dislocation diminishes accordingly. Additionally, an in- 
crease in temperature increases the dislocation mobility and speeds up the annihilation process, 
which results in an attendant increase in annihilation rates. The rate of annihilation is then modu- 
lated by the nucleation energy of a jog E^ og which can be calculated from atomistic simulations. 



3 Atomistic modeling of dislocations properties 

In the previous section we have identified the following set of material parameters required to 
estimate the contribution of each of the controlling unit processes: i? kink , L kink , U edge f/ screw 5 and 
Ei° g . In this section we briefly describe the computation of them using a First Principles-based 
Force Field with Molecular Dynamics. 

Quantum mechanics (QM) describes the atomic interactions from first principles, i.e. using no 
empirical input. Unfortunately QM methods are computationaly intensive and thus only applicable 
to small systems (hundreds of atoms) and short times (picoseconds). The study of most of the unit 
processes that govern the plasticity of materials (such as dislocation mobility, kink energies, etc.) 
involve many atoms and long simulation times. Such problems require the use of Force Fields, 
where the total energy of the system is given by a potential energy function of the atomic positions 
and does not involve the solution of Shrodinger's equation. The drawback of using potentials to 
describe the atomic interactions is that some accuracy is lost; it is thus of great importance to use 
accurate force fields to describe the atomic interactions. 



We developed a many body Force Field for Tantalum based on accurate QM calculations p4j ] 
which can be used with Molecular Dynamics (MD) to simulate systems containing millions of 
atoms. We fitted an Embedded Atom Model type Force Field (named qEAM FF) to a variety of ab 
initio calculations, including the zero temperature Equation of State (EOS) for bcc, fee, and A15 
phases of Ta in a wide pressure range, elastic constants, vacancy formation energy and energetics 
of a shear transformation in the twinning direction [p4|]. Ta is a bcc metal and no phase transition 
to other crustalline phase is known, but using QM we can calculate the EOS of thermodynamicaly 
unstable or metastable phases (such as A15, fee, hep, etc.). Including data about these high energy 
phases, with different coordination numbers, in the Force Field training set is important to correctly 
describe the atomic interactions near defects such as dislocations, grain boundaries, etc. 

We have used the qEAM with MD to study a variety of materials properties such as the melting 
temperature of Ta as a function of pressure [33p, dislocation properties [JT3J], and spall failure p5|]. 



In subsection [3J] we show the calculation of the core energy of edge and screw dislocations in 
Ta and in Subsection EO we calculate the double kink formation energy and nucleation length. 
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Figure 7: Differential displacement map of a relaxed quadrupole of screw dislocation in Ta. 

3.1 Core energy of the 1 /2a < 111 > screw and edge dislocations 

In order to study static properties of the l/2a< 111 > screw dislocation in Ta such us core struc- 
ture and energy we use a dislocation quadrupole in a simulation cell with periodic boundary con- 
ditions. Two of the dislocations have Burgers vector b=l/2a< 111 > and the other two have 
b=-l/2a< —1 — 1 — 1 >. Such an arrangement of dislocations minimizes the misfit of atoms on the 
periodic boundary due to the effects of periodic images. We build the dislocations using the atomic 
displacements obtained from elasticity theory and then we relax the atomic coordinates using the 
qEAM FF. In the bcc structure, there are two kinds of dislocation core configurations (easy core 
and hard core) that can be transformed to each other by reversing the Burgers vector [|TT|]. In this 
work we focus on the lower energy easy cores. In Figure [7] we show the differential displacement 
map (DDM) of our relaxed quadrupolar system. In the DDM atoms are represented by circles and 
projected on a (111) plane. The arrows represent the relative displacement in [111] direction of 
neighboring atoms due to the dislocation. We can see from Figure [7] that the equilibrium dislo- 
cation core obtained using qEAM FF has three-fold symmetry and spreads out in three < 112 > 
directions on 1 10 planes. 

Lets define strain energy as the total energy of our system once the perfect crystal energy is 
subtracted. The total strain energy can be divided in two terms: core energy (E c ) and elastic 
energy (E e ). The latter contains the self-energy of each dislocation and their interactions and can 
be calculated using linear elasticity theory. The core energy is the energy contained close to the 
dislocation line (closer than some distance r c called core radius) where, due to the large strain, 
elasticity theory is not valid and the details of the interatomic interactions are important. For our 
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Figure 8: Total strain energy of the quadrupolar system as a function of In d\/r c + A(di/d 2 ); the 
number of atoms in each simulation is shown. The dashed line is the linear fit to our atomistic data. 



quadrupole system the total strain energy takes the form [JO]]: 



E = E c (r c ) 



Kb 



In 



d\ 



A 



(14) 



where K depends on the elastic constants, di and d 2 are the nearest separation of dislocations along 
< 11 — 2 > and < 1 — 10 > directions and A(di/d 2 ) is a geometric factor which comes from the 
dislocation interactions. 

We studied quadrupolar dislocation cells of different sizes. In Figure (Q) we show the mini- 
mized energy as a function of \ndi/r c + A(d 1 /d 2 ) for the different simulation cells; in order to 
compare with previous calculations [JTTJ, [13]] we took r c = 2b. We can see that the total energies fol- 
low a straight line as predicted by elasticity theory (Eq. |14|), showing that the value chosen for the 
core radius is large enough to take account for the non elastic region near the dislocation line. From 
a linear fit to our data we determine the core energy E c = 1.30 eV/b and K = 3.33 x !0~ 2 eV/A3. 
The value of K can also be computed from the elastic constants giving 3.33 x lO~ 2 eV/ A3 in ex- 
cellent agreement with the one obtained from the fit. Recent ab initio calculations of core energy 
(using periodic cells containing 90 atoms) give 0.86 eV/b, lower than the value obtained with 
qEAM FF and the dislocation cores are compact and symmetric [JT3p . 

Using the qEAM we can calculate the strain energy associated with each atom. In Figure @ 
we show the atomic energy distribution (number of atoms per dislocation per Burgers vector as a 
function of their strain energy) for a system containing 5670 atoms in the periodic cell. We can 
see that there are 6 atoms with atomic strain energy higher than 0.15 eV and another 6 atoms with 
energy in the range 0.06-0.08 eV. They correspond to the 12 atoms per dislocation per Burgers 
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Figure 9: Histogram of atomistic strain energy distribution for the quadrupolar arrangement of 
screw dislocations. The cell contains 5670 atoms and is 7 Burgers vectors long. 



vector closer to the dislocation line and their total energy is 1 .35eV/b, very similar the core energy 
obtained from Eq. [14]. The rest of the atoms have lower strain energy and can be cosidered as the 
elastic part of the system. We can then define the dislocation core as formed by the 12 atoms per 
Burgers vector with higher energy. 

We have also calculated the core energy of the edge dislocation with 6=l/2a<lll>ona 
(110) plane. We build a simulation cell with axis oriented along < 112 > (x axis), < 110 > (y 
axis), and l/2a < 111 >(z axis); this cell contains 6 atoms. We then replicate the cell 3 times 
along x, 16 times along y, and 20 times along z; the number of atoms in the cell is then N=5760. 
We then remove 108 atoms to form a dipole of edge dislocations. Once the system is relaxed (both 
atoms and cell parameters) we have a 24.3967Ax75.1824.Ax56.632 Acell. Figure ( |TTj| ) shows a 
snapshot of the atoms projected on a < 112 > plane. 

In Figure [FT] we show the energy distribution for the edge dislocation (number of atoms per 
dislocation and per a < 112 > length as a function of their energy). Figure [IT] shows that the core 
of the edge dislocation contains atoms with higher energies and a broader distribution of energies 
as compared with the screw case [Figure @] . Taking into accound Figure [n] we define the core of 
the edge dislocation as formed by those atoms with strain energy higher that 0.1 eV. This definition 
leads to 36 atoms per a < 112 > or ~ 4.42 atoms per A and to a core energy of E e ^l = 0.827 
eV/A(in the case of the screw we had 12 atoms/b or ~ 4.17 atoms per A). The ratio between the 
core energy of the edge and that of the screw is: E^^ / E^l w ~ 1.77. It is important to mention 
that changing the number of atoms considered to belong to the core changes the core energy, but 
the difference is minor. Had we taken the 34 atoms per a < 112 > with higher energy as the core 
(leading to ~ 4.18 atoms /A, a density very similar to the one obtained in the screw dislocation) 
we get a very similar core energy: E e J^l = 0.80 eV/A. 
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3.2 Kink pair energy and nucleation lenght 



As already explained the kink pair mechanism controls the mobility of screw dislocations in bcc 
metals and atomistic simulations can provide the details of this mechanism. 

As we can see from Figure (|7|) the core of the screw dislocation spreads in three < 112 > direc- 
tions, this leads to two distint, but energetically equilavent, core configurations; we will name them 
positive (P) and negative (N) cores p6|]. The shorter (and lower energy) kinks possible involves the 
displacement of the position of the dislocation line in the (111) plane from one equilibrium posi- 
tion to a nearest neighbor equilibrium position; the displacement involved isl/3a<112>. There 
are six possible < 112 > directions but only two need to be considered by symmetry, this leads to 
two kink directions which we will call left (L) and right (R). The two dislocation cores (N and P) 
and two directions (L and R) lead to 8 different single kinks: NRP, NRN, PRP, PRN, NLP, NLN, 



PLP and PLN. We have studied all of them in detail [|46|], here we will concentrate on the single 
kinks that lead to the lowest energy kink pair. We calculated the formation energy and length of the 
various kinks using quadrupolar arrangements of dislocations as explained in subsection BTT1. The 



simulation cell lengths are 40.7 A in the [11-2] direction, 42.3 A in [1-10] and 431.8 A in [111] 
containing 40500 atoms; the details of these calculations can be found in PBfl, We calculate the 
kink energy as the difference of strain energy between the quadrupolar systems containing kinks 
and that corresponding to perfect straight dislocations. This energy difference is divided by four 
to get the energy per kink. Using the qEAM FF we find that the lowest energy kink pair is formed 
combining the PLN and NRP kinks. We define the kink pair nucleation energy as the sum of the 
formation energy of the two songle kinks leading to E k%nk = 0.725 eV. The nucleation energy 
calculated in this way does not take into account the attractive interaction between the two kinks 
which lowerers the nucleation energy. This interation energy is very small for separation larger 
than~ 15 b IJH^j. 

As explained above, a critical parameter for the micro-mechanical modeling of plasticity is, 
apart from the kink pair energy, its nucleation length L kink . We studied both the energetics and 
structure of the variuos kinks along the dislocation line. Figure [12] shows the extent of the kinks 
both from structural and energetic points of view. We show the position of the dislocation in the 
direction of the kink along the dislocation line for a PLN kink [Figure ]I2](a)] and NRP kink [Figure 
|T^|(c)]. We also show the total energy of the quadrupolar system along the dislocation line for the 
PLN [Figure |i~2|(b)] and NRP [Figure |i~2|(d)] kinks. This is calculated by dividing the system in 
the [111] direction in regions of width equal to the Burgers vector and calculating the total energy 
in each slice. The structural length of the PLN kinks is L%£ N = 8 b [Figure [12] (a)]; while its 
"energetic extent" is L^ e N = 14 b [Figure [T2| (b)] . For NRP kinks we obtain: Lf t ^ p = 8 b [Figure 
g (c)] and L^ p = 20 b [Figure p| (d)] . 



Going back to the definitions of the paramenters entering the equation that governs the dislo- 
cation mobility [Equations ([[]) and @]; the effective Peierls stress (r ), Equation (Q), is defined as 
the applied stress for which the nucleation free energy for a kink pair (AG) is zero. AG is given 
by: 



AG = E kink ± Tl P bL kink , (15) 

where L kmk is the effective kink pair nucleation length and lp is distance advanced by the disloca- 
tions; in the kinks studied here lp = \l/3a < 112 > |. The second term in the right hand side of 
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Figure 12: PLN and NRP kinks in Ta using the qEAM FF. (Top Left) PLN kink: Dislocation 
position in the [11-2] direction along the dislocation line; we can see the dislocation moves from 
one equilibrium position to the next in a distance of 10 Burgers vectors. (Top Right) PLN kink: 
total energy in the quadrupolar system with four PLN kinks along the dislocation line. The system 
is divided in slices with thickness equal to b and the energy in each region is calculated. (Bottom 
Left) NRP kink: Dislocation position in the [11-2] direction along the dislocation line; we can see 
the dislocation moves from one equilibrium position to the next in a distance of 10 Burgers vectors. 
(Bottom Right) NRP kink: total energy in the quadrupolar system with four PLN kinks along the 
dislocation line. 
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Figure 13: Schematic diagram of a kink pair formed by a NRP and PLN single kinks. The four 
terms entering in the work expresion (Equation [TJ) are shown in the Figure. 



Equation [15] is the work done by the external stress when the kink is nucleated. Figure [13] shows a 
schematic diagram of a PLN-NRP kink pair. We can see that the work done by the external stress 
to nucleate the kink pair can be divided in four terms: 

/ r PLN j PLN _ j PLN j NRP _ j NRP j NRP \ 

rbl P L kink = rblp + s -^— + + (16) 



where L n is the effective kink pair length. In Figure IT3] we show the four terms in the right hand 



side of Equation (|16|). Note that Equation [16| assumes that the kinks are straight lines connecting 
the two equilibrium positions of the dislocation. In this way we obtain the effective kink pair 
nucleation length L kink = 17 b. 

The remaining materials parameter is the nucleation energy of a jog E*° s . In this work we take 
E* 0B as the PLN-NRP kink pair nucleation energy. 



4 Experiment, Validation and Prediction 

To test the predictive capabilities of the multiscale approach we first select a set of material parame- 
ters to best fit the experimental results, then we compare these parameters against the atomistically 
computed ones and finally we predict the macroscopic response using the atomistics parameters. 
As we shall see, the agreement between the fitted and computed by atomistics material parameters 
is remarkable, and the macroscopic predicted response retains most of the experimental features. 
These facts provide confidence in the multiscale modeling approach, indicating that even in the 
case that experimental data would not have been available, still the macroscopic behavior could 
have been predicted based only on atomistic calculations. 

The experiment data correspond to uniaxial tests on Ta single crystals of Mitchell and Spitzig 
[p7|]. In these tests, 99.97%-pure Ta specimens were loaded in tension along the [213] crystallo- 
graphic axis, at various combinations of temperature and strain rate. In particular we considered 
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Table 2: Material parameters for Tantalum 
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temperatures ranging from 296 K to 573 K, and strain rates ranging from 10 _1 s~ l to 10~ 5 s -1 . 
The numerical procedure employed for the integration of the constitutive equations has been de- 
scribed elsewhere [^8]]. The constitutive update is fully implicit, with the active systems determined 
iteratively so as to minimize an incremental work function. All stress-strain curves are reported in 
terms of nominal stress and engineering strain. 

Two different sets of material properties were used for the numerical simulations. The first set 
was obtained by fitting the simulation results to the experimental results. Table || identifies the sub- 
set of parameters which are also amenable to direct calculation by atomistic based methods. The 
table lists the parameter values obtained by these methods, as described in the previous sections, 
in parallel with the values obtained by the fitting approach. Thus, in the second set of properties 
which was used for numerical simulations, atomistic-based values replace fit-based values, when 
available. This is the case for the edge and screw dislocation self-energies, as well as the kink-pair 
formation energy and length. Clearly, those two sets don't differ by much, which strongly support 
the validity of the advertised multiscale paradigm. For a complete list of parameters for the model, 
the reader should refer to [|7p . 

Figs. [R| and [T3| show the predicted and measured stress-strain curves for a [213] Ta crystal over 
a range of temperatures and strain rates. One can compare, from top to bottom: the experimental 
results, the results obtained after fitting the parameters, and the results obtained with atomistic- 
based parameters. It is evident from these figures that the model, with both sets of parameters, 
captures salient features of the behavior of Ta crystals such as: the dependence of the initial yield 
point on temperature and strain rate; the presence of a marked stage I of easy glide, specially at 
low temperature and high strain rates; the sharp onset of stage II hardening and its tendency to 
shift towards lower strains, and eventually disappear, as the temperature increases or the strain rate 
increases; the parabolic stage II hardening at low strain rates or high temperatures; the stage II 
softening at high strain rates or low temperatures; the trend towards saturation at high strains; and 
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(b) Predictions of the model with fitted parameters 
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(c) Predictions of the model with atomistic parameters 

Figure 14: Temperature dependence of stress-strain curves for [213] Ta single crystal (e 
s _1 ). 20 
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(b) Predictions of the model with fitted parameters 
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(c) Predictions of the model with atomistic parameters 

Figure 15: Strain-rate dependence of stress-strain curves for [213] Ta single crystal (T=373 K). 

21 



the temperature and strain-rate dependence of the saturation stress. Thus, the predictive approach 
based on atomistic methods clearly shows its capacity to produce results matching the experimental 
evidence. 

The theory reveals useful insights into the mechanisms underlying these behaviors. For in- 
stance, since during state I the crystal deforms in single slip and the secondary dislocation densities 
are low, the Peierls resistance dominates and the temperature and strain-rate dependency of yield 
owes mainly to the thermally activated formation of kinks and crossing of forest dislocations. It is 
interesting to note that during this stage the effect of increasing (decreasing) temperature is similar 
to the effect of decreasing (increasing) strain rate, as noted by Tang et al. [|2"3]]. The onset of stage 
II is due to the activation of secondary systems. The rate at which these secondary systems harden 
during stage I depends on the rate of dislocation multiplication in the primary system. This rate is 
in turn sensitive to the saturation strain 7 sat , which increases with strain rate and decreases with 
temperature. As a result, the length of the stage I of hardening is predicted to increase with strain 
rate and decrease with temperature, as observed experimentally. Finally, the saturation stress is 
mainly governed by the forest hardening mechanism and, in particular, by the strength of the forest 
obstacles. This process is less thermally activated than the Peierls stress, since the corresponding 
energy barriers are comparatively higher. Consequently, the stress- strain curves tend to converge 
in this regime, in keeping with observation. 

The apparent softening observed in simulation results at the lowest temperature (296 K) and 
the highest strain rate (1CT 1 s _1 ) is actually an effect of the boundary conditions, allowing some 
level of rotation of the specimen. Since in those cases, the material hardening is relatively low 
(stage I only), this geometrical softening dominates in the apparent macroscopic behavior. In the 
other cases, the activation of several systems at high strains results in a more isotropic deformation, 
in turn leading to limited rotations. In order to take the exact experimental boundary conditions 
into account, a finite element model of the whole specimen should be used, allowing for a non- 
homogeneous deformation field. 
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